Scaling of avalanche queues in directed dissipative sandpiles 



Bosiljka Tadic 1 '* and Vyatcheslav Priezzhev 2 '** 
1 Jozef Stefan Institute, P.O. Box 3000, 1001 Ljubljana, Slovenia 
2 Bogolubov Laboratory of Theoretical Physics, Joint Institute of Nuclear Research, 14-1980 Dubna, Russia 



Using numerical simulations and analytical methods we study a two-dimensional directed sandpile 
automaton with nonconservative random defects (concentration c) and varying driving rate r. The 
automaton is driven only at the top row and driving rate is measured by the number of added 
particles per time step of avalanche evolution. The probability distribution of duration of elementary 
avalanches at zero driving rate is exactly given by Pi(t, c) = i~ 3//2 exp [t ln(l — c)]. For driving rates 
in the interval < r < 1 the avalanches are queuing one after another, making increase the periods 
of non-interrupted activity of the automaton. Recognizing the probability Pi as a distribution of 
service time of jobs arriving at a server with frequency r, the model represents a new example of the 
(E, 1, GI/oo/1) server queue in the queue theory. We study scaling properties of the busy period 
and dissipated energy of sequences of non-interrupted activity. In the limit c — > and varying linear 
system size L <C 1/c we find that at driving rates r < L -1 ^ 2 the distributions of duration and 
energy of the avalanche queues are characterized by a multifractal scaling and we determine the 
corresponding spectral functions. For i>l/c increasing of the driving rate somewhat compensates 
the energy losses at defects above the line r ~ y^c. The scaling exponents of the distributions in 
this region of phase diagram vary approximately linearly with the driving rate. Using properties of 
recurrent states and the probability theory we determine analytically the exact upper bound of the 
probability distribution of busy periods. In the case of conservative dynamics c = the probability 
of a continuous flow increases as F(oo) ~ r 2 for small driving rates. 



I. INTRODUCTION 

In the past decade the sandpile type of cellular au- 
tomata played a special role in understanding of self- 
organized criticality in nonlinear dynamical systems (for 
a recent review see Ref. Q). In sandpile automata the 
properties of the dynamics which are essential for the 
occurrence of self-organized critical states can be moni- 
tored in a direct manner. Apart from the relaxation rules, 
these are the following properties: type of driving and 
time-scale separation; conservation law of the dynamics; 
direction of mass flow; role of boundaries. In addition, 
the Abelian nature of the toppling rules in some sandpile 
automata enabled derivation of certain exact results , 
in contrast to other dynamical systems where such cal- 
culations are not available. Numerous sandpile models, 
both deterministic and stochastic Jj| are shown to exhibit 
dynamic critical states in the limit of "infinitely slow" 
driving (i.e., at zero driving rate r = 0). In this limit 
a new avalanche is initiated only after the previous one 
has stopped, thus the time scale separation is exactly ob- 
served. On the time scale of perturbations, the avalanche 
evolution is seen as occurring instantly. The existence of 
the critical states in the case of directed Abelian sandpile 
automaton at zero driving rate has been proved exactly 
by Dhar and Ramaswamy . At this point it is interest- 
ing to mention that the model studied in Ref. Q is char- 
acterized by local driving and deterministic conservative 
dynamics. The automaton with conservative stochastic 



dynamics, on the other hand, has been shown to belong 
to another universality class Q . The presence of noncon- 
servative defects in Dhar-Ramaswamy automaton leads 
to a subcritical behavior (jjj . 

The behavior of driven dynamical systems at finite 
driving rates (r > 0) represents an important subject 
both for theoretical and practical reasons. A finite driv- 
ing rate may appear either as a control parameter set 
from outside, or as a probability distribution originating 
from another coupled stochastic process. In practice, the 
systems are driven by an external field, which oscillates 
with a finite frequency. Examples are Barkhausen noise 
Q , integrate and fire oscillators |?J , granular material in 
rotating drums etc. Queuing jobs at a server [||, 
e.g., in teletraffic, is an example where the frequency of 
arriving jobs is given by a random process. 

Theoretically at finite driving rates r > the probabil- 
ity that a new avalanche starts before previous one has 
stopped increases with increasing r. This obviously leads 
to different statistics of avalanches, where an avalanche 
is understood to represent a non-interrupted activity of 
the system. For large driving rates a continuous flow (an 
avalanche which never stops) is expected in sandpiles. 
Similarly, a single spanning cluster may occur in driven 
disordered systems. Therefore, a time scale separation 
becomes less and less apparent with increasing r. In ad- 
dition, by increasing driving rates, the local driving looses 
its strict sense. Thus fast driven sandpiles are placed be- 
tween strictly local driving, where the system is driven 
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at a single (random) site, and global driving, where the 
same perturbation applies to all sites in the system. The 
role of the conservation law (conservation of number of 
grains in the interior of the system) is also expected to 
be changed at finite driving rates. In r — > limit, locally 
driven nonconservative systems appear to be subcritical 
JlO| , whereas when the driving is global the critical states 
may appear even if the dynamics is dissipativc ]ll|,f^| . So 
far neither an exact theory nor a renormalization-group 
analysis of fast driven critical systems has been done. A 
general questions as the existence of critical states at fi- 
nite driving rates and universal scaling properties of the 
system in the limit of large distances and long times re- 
main yet to be understood. 

Recently two numerical simulations elucidated certain 
important properties of sandpiles at finite driving rates. 
In a 1-dimensional ricepile model Corral and Paczuski 
fiJI have first introduced avalanches for a finite driv- 
ing rate as non-interrupted periods of activity and have 
shown that these avalanches diverge at rates r > r c (L) ~ 
i" 20 for a given L. The rational behind this conclusion 
is that an ever-running avalanche occurs for driving rates 
r ~ l/(t) , where (f) ~ L z{2 ~ r ^ is the average duration 
of avalanches in zero driving rate ]l2|] . In another exam- 
ple Barrat et al. [ ft3| have shown that in a 2-dimensional 
symmetric Abelian sandpile model mixing of time scales 
at finite driving rates leads to correlations which appear 
to violate the fluctuation-dissipation theorem. 

In this work we study a simple 2-dimensional model 
with strictly directed flow of grains and deterministic 
toppling rules. We add particles only at the top row 
with driving rate r. The driving rate r is defined as 
a number of added particles per time step of avalanche 
evolution. We consider both conservative and nonconser- 
vative dynamics. A fraction of sites c are considered as 
annealed nonconservative defects. By toppling at a defect 
site two grains are lost, thus affecting the propagation of 
avalanche below that site 0] . When c = the dynamics 
is strictly conservative. In the r = limit and c = the 
model has been exactly solved by Dhar and Ramaswamy 
H . The critical states where shown to consist of heights 
h = and h = 1 occurring with equal probability 1/2 at 
each lattice site. The duration of an avalanche is given 
by the probability distribution for large t 



Pi(t) ~* 



-3/2 



(1) 



In the presence of nonconservative defects it has been 
shown [p|JT^] that the screening of the power-law distri- 
bution in Eq. (fy) occurs as 



P 1 (t)~t- 3 / 2 exp(-t/0 ; £~l/c 



(2) 



In this paper we will refer to the distributions in 
Eqs. (0) and (^) as probability distributions of elemen- 
tary avalanches, to be distinguished from the combined 
avalanches, which occur at finite driving rates and which 



consist of series of elementary avalanches. Up to rela- 
tively high driving rates r = 1 the model has the prop- 
erty that successive elementary avalanches are running 
one after the other, in contrast to cases studied in Refs. 
p^LpI, where merging of avalanches may occur at any 
finite r. After an elementary avalanche is over the sys- 
tem is characterized by statistically unchanged distribu- 
tion of heights, owing to a weak correlation between the 
avalanches in the recurrent states. A finite probability 
of avalanche collision, which accelerates flow of grains, 
occurs in this model only for r > 1. Here we restrict the 
study to the case r < 1 , where the formation of avalanche 
queues is a dominant phenomenon which determines the 
scaling properties of the system. 

The problem of avalanche queue in our model can be 
regarded as an example of the server queue of the 
class (E, 1, GI/oo, 1), which is studied as a model in the 
analysis of various applied problems, e.g., in telecommu- 
nications, insurance, etc. These analogies are made clear 
by noticing that the avalanches of a directed sandpile 
model have the corresponding terms in the language of 
the queue theory as follows: (i) elementary avalanche — 
customer; (ii) duration of elementary avalanche — ser- 
vice time; (iii) driving rate — frequency of arrivals of 
customers; (iv) number of elementary avalanches coex- 
isting at a given moment of time — number of working 
servers; (v) duration of a combined avalanche — busy 
period. The notation (E, 1, G//oo, 1) means that we 
deal with customers arriving by one and being served 
by one. The letter E means that the arrival times are 
generated by a Bernoulli process with the distribution 
Prob(t = k) = p k {l-p);k = 0,l,2,...;p > . The 
symbol GI/oo means that the service times are identical 
independently distributed (i.i.d.) random variables and 
the infinite number of servers provides a non-restricted 
number of customers which can be served simultaneously. 
Despite a huge literature devoted to this subject ||, most 
of papers focus on the distribution of q n , the number of 
working servers at a given moment of time. The impor- 
tance of the tail behavior of the busy period distribu- 
tion for fluid queues in telecommunications, generalized 
processor sharing and other applications was recently 
pointed out in Ref. |lq| . Here we concentrate on some 
other properties of the queue: the scaling behavior of the 
busy period and dissipated energy distribution. 

Given the duration of elementary avalanches in Eq. (p]) , 
we may conclude that the directed sandpile model for 
r < 1 represents a special case of the queue theory with 
the power-law distribution of the service time with the 
exponent v = 3/2 — 1 = 1/2. This implies that the aver- 
age service time per customer diverges (t) — ► oo. In prac- 
tice, service times are restricted to finite values, which 
corresponds to the power-law distribution with 1 < v < 2 
p6[ . This may explain why the queue with the distribu- 
tion of the type given in Eq. (jlj) has not been studied so 
far. In our model a finite average duration of elementary 
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avalanches (t) < oo is achieved in two cases: (a) In the 
case of distribution in Eq. ([!]) when the system size L is 
finite, hence the distribution is truncated at t = L; (b) 
In the case of finite dissipation c > 0, where the distri- 
bution in Eq. (^) has a characteristic duration £ < oo for 
all finite c values (see Refs. [||,|l5| and below). 

For a finite L we find a continuous flow phase (F) for 
low dissipation and large driving rates, and three regions 
with intermittent behavior of avalanche queues. These 
are regions with subcritical (S), nonuniversal (N), and 
multifractal (M) behavior, shown schematically on the 
phase diagram in Fig. 1. When the length separation 
L> ( holds, we find a line in the (r, c)-plane where loss 
of particles on defects in the interior of the pile becomes 
"compensated" by fast adding of particles from outside. 
In the region above the compensation line the avalanche 
queues exhibit a scaling behavior with the scaling expo- 
nents depending on the driving rate: The slopes decrease 
whereas the fractal dimensions increase with driving rate. 
Cut-offs with a stretch-exponential behavior appear. In 
the limit c — > and when the system size L -C £ is varied 
multifractal scaling properties describe the scaled distri- 
butions, rather than a simple finite-size scaling. 




Concentration of Defects 

FIG. 1. Schematic phase diagram for a finite system size 
L. Cross-hatched area represents a crossover region between 
(right) and £ 3> L (left). Regions with distinct behav- 
ior of the avalanche queues are shown: subcritical (S), nonuni- 
versal (N), and multifractal scaling region (M), and flow phase 
(F). Solid line represents the compensation line. Transition 
to the flow phase is marked by o for c = and by dotted line 
for small c > 0. In the origin only single Dhar-Ramaswamy 
avalanches occur. 

In general, the distributions of combined avalanches 
are characterized by a scaling function of two arguments 
of the form 

P(X, r, L c ) ~ X-^T(XL- D * , rL l J 2 ) , (3) 

where X represents either duration, t, or number of top- 
plings fT^ ], n, and L c = min (£, L). The corresponding 
fractal dimensions Dx are defined by 

(x)t ~ e Dx , (4) 



where the average is taken over all combined avalanches 
of a selected length I measured along the direction of 
transport. 

Using analogy to the queue theory and the properties 
of the recurrent states we were able to derive an exact 
upper limit of the distributions of busy periods and to 
discuss the limit L — > oo. We also derived the expression 
for the probability of continuous flow in the conservative 
limit. 

The paper is organized as follows: In Sec. II we define 
the model and consider the case of conservative dynamics 
by numerical simulations on finite lattice. In Sec. Ill we 
present results of simulations in the case of finite concen- 
tration of nonconservative defects. In Sec. IV we present 
details of the analytical results. The paper contains a 
short summary of the results and discussion in Sec. V. 



II. MULTIFRACTAL QUEUES OF 
DHAR-RAMASWAMY AVALANCHES 

The sandpile automaton model introduced by Dhar 
and Ramaswamy represents an example of a self- 
organized criticality with exact solution in the limit of 
zero driving rate ^ . In this Section we consider the same 
model at finite driving rates < r < 1. The dynamic 
rules of the model are summarized as follows : We con- 
sider a 2-dimensional square lattice oriented downwards, 
with a dynamic variable — height h(i,j) — associated at 
each site. Grains are added at the top row only, and 
mass flow is only down. The toppling at a site oc- 
curs deterministically whenever h(i,j) ~>h c = 2, and two 
grains are transferred downward, i.e., 

h(i,j)^h(i,j)-2; &(* + l,j±)->ft(* + l,j±) + l , 

(5) 

where (i + l,j±) represents two downward neighboring 
sites to the site 

The probability distribution of duration of avalanches 
in zero driving rate -Pi(t) ~ t~ Tt V(tL~ v ) with r t ° = 3/2 
given in Eq. (|l|) becomes exact at large t . In addition 
the dynamic exponent 2 =1. This implies that the av- 
erage duration (t)o ~ L 1 / 2 — > oo. Similarly, the area s 
enclosed by the boundary of an avalanche is given by the 
distribution D(s,L) ~ s- r °V(sL- D °), where r s ° = 4/3 
and the fractal dimension D° s =3/2. Note that the num- 
ber of toppled grains at each active site is two, then the 
distribution of the number of toppled grains within an 
avalanche, D(n), is described by the same exponents, i.e., 
r° =4/3 and D° = 3/2 at zero driving rate. 

A finite driving rate r > is implemented as follows. 
An avalanche is started from the top and at each step of 
the avalanche progress a new particle is added with prob- 
ability rata random site at first row. We also consider 
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a deterministic addition of particles, i.e., we add a parti- 
cle at regular intervals At. Both approaches lead to the 
same results when the statistics is high enough. An added 
particle may trigger a new elementary avalanche before 
the previous one stops, thus making a pattern of active 
sites distributed over lattice. A snapshot of growth of a 
combined avalanche with marked active sites is shown in 
Fig. 2 (top) . A combined avalanche — avalanche queue — 
is thus determined by a non-interrupted activity on the 
lattice and it stops when no more active sites occur. Then 
a new avalanche is started. It should be noted that when 
r > the number of added grains is higher than the 
number of combined avalanches. Another important re- 
mark is that the repeated toppling at a site may occur 
as soon as r > 0. This leads to the inequality (n) > (s), 
and thus D n > D Sl and z > 1. The numerical simula- 
tions confirm these conclusions (see below). Typically, 
we consider 2 x 10 6 combined avalanches at each driving 
rate and lattice size. Periodic boundary conditions are 
applied in the perpendicular direction. 




FIG. 2. Top: Growth of a combined avalanche in the case 
of conservative dynamics (c = 0) for driving rate r — 0.05 
and L — 100. Nine fronts of active sites (dark) are visible. 
Bottom: A complete combined avalanche in the case of dis- 
sipative dynamics with c = 0.02 and driving rate r = 0.5. 
Different colors correspond to distinct toppling waves. 

In this Section we perform numerical simulations for 
r > and finite lattice sizes L. The limit L — > oo 
will be discussed in Sec. IV. In main Figs. 3 and 4 are 
shown the integrated probability distributions of dura- 



tion P(t! > t,L) and number of topplings D(n' > n, L) 
for fixed driving rate r = 0.05 and various lattice sizes L. 
It should be noted that both slopes and cut-offs of these 
distributions appear to be different compared to ones of 
the elementary avalanches. In particular, slopes decrease 
with r (see more detailed discussion in next Section). For 
instance at r — 0.05 we find r t = 0.4 and r n = 0.3, in the 
steep part, and r t = 0.31 and r n = 0.2 in the flat part 
near the cut-off. A cut-off in the probability distribution 
of durations appear (cf. main Fig. 3). The characteristic 
jump at t = L is related to the conditional probability: 
an activity lasts longer than L steps only if the preced- 
ing avalanches last not shorter than t = L. The jump 
decreases with increasing lattice size. 
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FIG. 3. Double logarithmic plot of the integrated probabil- 
ity distribution of duration (busy periods) P(t, L) of queues 
of Dhar-Ramaswamy avalanches vs. duration t, measured in 
Monte Carlo steps [MCS]. Fixed driving rate r = 0.05 and 
various lattice sizes L =24,48,96,192, and 384 are used. Inset: 
Multifractal spectral function $ t (at) vs. at obtained from 
the data in main figure according to Eq. (6) using Xo = 1, 
L = 2. 

It is interesting that these distributions can not be 
scaled according to a simple finite-size scaling (with new 
exponents), as one may naively expect. Instead, we find 
that a multifractal scaling applies according to the law 



P(X,L,r)=[- 



*xOx) 



a x 



log(X/X a ] 
log(i/i ) 



(6) 



where, as before, X stands for t or n. The corresponding 
spectral functions $ t (cv t ) and <E>„(a„) are determined nu- 
merically for r = 0.05 and shown in the insets to Figs. 3 
and 4, respectively. The spectrum depends on the driv- 
ing rate. For driving rates close to the line r ~ L^ 1 / 2 
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extremely large avalanches may appear and the scaling 
fits fail. The origin of multifractality in the queues of 
Dhar-Ramaswamy avalanches can be found in the fact 
that an unrestricted multiple toppling at each site may 
occur, and that a toppling at a given site releases a local 
avalanche which propagates from that site downwards. 
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FIG. 4. Same as Fig. 3 but for the integrated distribution 
of number of topplings (mass) D(n, L) vs. mass n [number 
of grains]. Shown are only curves for L =96,192, and 384. 
Inset: Spectral function $„(q„) vs. q„. Xq = 0.97 ± 0.02, 
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2,15 ±0.05. 
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FIG. 5. Attempted finite-size scaling fit of the data from 

main Fig. 3 and 4 for the distribution of avalanche duration 
(right) and mass (left). 

For comparison we show how a finite-size scaling fit of 
the data fails. In Fig. 5 we present an attempt of scaling 
collapse of the data shown in Figs. 3 and 4 above, ac- 
cording to the formula P(X, L) = L- X V{X L~ D x). Note 
that the best collapse of the tails of distributions are ob- 
tained by fixing the fractal dimensions as D t — 1.20 and 
D n = 1.75, which correspond to ax at the shoulder of 
the spectrum $t(at) and $„(a„), respectively. Fixing 



a smaller (larger) value for the fractal dimension leads 
to systematic shifts of the distribution tails to the right 
(left) with increasing L. The best fit shown in Fig. 5 is 
obtained for A = 0.35, which satisfies (within numerical 
accuracy) the scaling relation A = z(rt — 1) = D n (T n — 1) 
with tx — 1 determined at the flat part of the distribu- 
tion. However, as the Fig. 5 shows, fixing Dx and A 
leads to the systematic shifts of the 'horizontal' part of 
the distributions to the right with decreasing L. Fixing 
the exponents independently from the scaling relation re- 
sults in crossing of the lines for different L values. 

For driving rates r > L -1 / 2 an ever-running avalanche 
may occur, representing a continuous activity on the lat- 
tice. The flow phase can be characterized by an average 
number of topplings per site, which is expected to have a 
nontrivial L-dependence. The probability of occurrence 
of the flow phase in the limit L — > oo will be discussed in 
Sect. IV. 



III. NONUNIVERSALITY IN DISSIPATIVE 
DYNAMICS 

In the presence of dissipative defects c > the distri- 
bution of elementary avalanches, which is given in Eq. 
(Q), appears to have a finite characteristic length £ < oo. 
Thus the average duration at zero driving rate is finite 
(i)o < oo. Precise value of the average duration is con- 
trolled by an external parameter — probability of dissipa- 
tion c, and not by system size L, provided that L > £. 
The screening length £ ~ 1/c was first estimated numeri- 
cally in Ref . Q . A more precise analytical expression can 
be derived (see below and Ref. fllslO as ~ — ln(l — c). 
Here we perform numerical simulations in the case L > £ 
at driving rates < r < 1. In this range of driving rates 
we expect the role of lattice size in the analysis of Sec. 
II is to be replaced by the characteristic length £. In ad- 
dition, competition between dissipation and driving rate 
leads to new phenomena. 

In Fig. 2 (bottom) an example of a combined avalanche 
is shown for c = 0.02 and driving rate r = 0.5. It is re- 
markable that the number of topplings per site decrease 
with distance from the top. Intermittency of the dynam- 
ics as well as the occurrence of the long-range correlations 
can be seen by direct examination of the recorded activ- 
ity of the system n(t) at each time step. For the server 
queue the quantity n(t) is interesting as the measure of 
the energy which is dissipated by the server at a given 
moment of time t. In Fig. 6 we show an example of the 
recorded signal for certain choice of parameters r, c, and 
L corresponding to the region (N) of the phase diagram 
(cf. Fig. 1). A combined avalanche on this recording 
is represented by a set of peaks between two consecu- 
tive drops of the signal to the base line. The Fourier 
spectrum of the signal (shown on top panel in Fig. 6) 
exhibits a power-law behavior. The slope if ~ 0.9 weakly 
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increases with driving rate r. 
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FIG. 6. Sandpile noise n(t) — number of active sites at 
time t, plotted against time t [MCS] (lower panel) and its 
Fourier spectrum (top panel) for fixed driving rate r=0.33. 
Dissipation rate c = 0.01 and L — 128. Shaded area shown 
in the lower panel is an example of unit signal, corresponding 
to a combined avalanche. 



of values of driving rates the scaling exponents decrease 
and the fractal dimensions increase with r, while the scal- 
ing relations between various exponents are found to be 
satisfied within numerical error bars. 
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FIG. 7. Distribution of avalanche mass D(n, r) vs. mass n 
[number of grains] shown in double logarithmic scale for fixed 
dissipation rate c = 0.01 and L — 192 and for various driving 
rates r= 0.021,0.041,0.083,0.125,0.166,0.25, and 0.33 (left to 
right). Inset: Average mass (top), area (middle) and duration 
(bottom) of combined avalanches plotted against driving rate 
[MCS -1 ] for fixed c = 0.01 and L — 192. At each point aver- 
age is taken over 2 x 10 6 combined avalanches. Fitting curves: 
(n) = 167 x exp(11.4r 128 ), and (t) = 23 x exp (9.6r 142 ). 



The scaling properties of the distributions of avalanche 
queues depend on the mutual ratio of the driving and dis- 
sipation rates. In particular, the scaling function in Eq. 
(||) exhibits a nontrivial dependence on both arguments 
x = and y = r£}/ 2 . (Another suitable choice of 

variables would be (tc, rt 1 / 2 ).) In Fig. 7 we show the dis- 
tribution of the avalanche mass (dissipated energy) n for 
fixed c = 0.01 and various driving rates r. In general, the 
cut-offs of the distributions increase and slopes decrease 
with increased driving rates r. More detailed analysis of 
the slopes shows that for the range of values of driving 
rates the scaling behavior of the distributions can be de- 
scribed by the scaling exponents which depend on the 
driving rate. The slopes of various distributions and the 
corresponding fractal dimensions, which are defined in 
Eq. (0), are shown against driving rate in Fig. 8. The 
dissipation rate is fixed c = 0.01. Note that T£ in Fig. 
8 represents slope of the distribution of the largest lin- 
ear length reached in a combined avalanche. For a range 



The variation of the scaling exponents can be approx- 
imated by a linear dependence of r. It is interesting to 
note that a qualitatively similar behavior — linear varia- 
tion of the scaling exponents with driving rate, has been 
measured experimentally in the case of Barkhausen noise 
in driven disordered ferromagnets Q. Our present anal- 
ysis suggests that such behavior can be related to an 
interplay of driving rate and dissipation at defects, and 
that it applies more universally. Dependence of the frac- 
tal dimensions (and of the slope exponents via scaling 
relations) of the avalanche queues on driving rate r can 
be linked to the r-dependence of the average length of 
the queue < N > = 1 + r < t >q. It has been shown 
recently [jOI that the fractal spectrum of the series of 
elementary signals in the case of transit times in the ri- 
cepile model varies as a power of the length of series. 
A precise r-dependence of the exponents in the case of 
avalanche queues requires additional work and will be 
given elsewhere. 
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FIG. 8. Lower panel: Various scaling exponents of the 
avalanche queues plotted against driving rate for fixed dis- 
sipation c = 0.01 and L = 192. * indicates the products 
z(r t - 1) « D„(t„ - 1) « D 3 (r s - 1) « l.l(r/ - 1). Also 
shown are linear fits of the data. Top panel: Amplitudes 
£n (upper curve) and £t defined in Eq. (7) vs driving rate. 
Lines: fit curves satisfying Eq. (8) with B n — 6.29 ±0.46 and 
B t =4.5 ±0.51. 

The opposite effects on the exponents are obtained by 
increasing concentration of defects c at fixed driving rate. 
In particular, the slopes of the distributions increase and 
fractal dimensions decrease with increasing concentration 
of defects in a limited range jlj| . 

An exact expression for the scaling function V(x, y) in 
Eq. (||) can not be guessed. It appears that the cut of 
the surface V{x 1 y) at r — const for well balanced values 
of c and r can be approximated by a stretch-exponential 
function, so that we have 

P(X, C) r) = X-^Mexp(-X CT */6Kr)) ■ (7) 

Here X stands for n or i, and we find a n = 1.14 ± 0.04 
and at = 1.28 ± 0.06, for the distribution of energy and 
duration, respectively. The amplitudes £jr (v) can be fit- 
ted (see top panel of Fig. 8) by the following function of 
driving rate r 

Zx(r)=A x (c)exp(rB x (c)), (8) 

for a fixed dissipation rate c. 

The observed parameter-dependence of the probability 
distributions is also reflected in the behavior of the av- 
erage duration and energy of combined avalanches. No- 
tice that the average duration (t) represents the average 



busy period of a server in the queue theory. Beside the 
average duration (t) , we also compute the average values 
of number of topplings (energy) (n) and area (s) (total 
number of distinct sites) affected by the processing of a 
combined avalanche. These average values are shown vs. 
r in the inset to Fig. 7 for c = 0.01. The values (t) and 
(n) increase with driving rate faster than an exponential 
function. Fitting the data in the inset to Fig. 7 we find 

(X) = a xcxp(aix(c)r CT ) , (9) 

where X = n,t and we estimate a = 1.2 ± 0.1. Note 
that, in contrast to durations and energies, the average 
area of an avalanche is bounded by the number of cells 
in the system (s) < L 2 . 

As already mentioned above, for a given dissipation 
rate c and L»(a driving rate r^ic) exists such that fast 
addition of grains compensates the losses in the bulk. In 
fact along an extremal line 

ro(c) ~ Kf 1/2 » kVc (10) 

the coherence length remains constant. The existence of 
the compensation line Eq. ( |To| ) can be demonstrated by 
considering sets of data for average durations and ener- 
gies against driving rate r, obtained for different charac- 
teristic length £ ~ 1/c. These data can be scaled accord- 
ing to the following scaling form 

( X ^-D x (2-rx) = g ^*(2-n) _ ^ _ (11) 

The corresponding scaling fits for the two cases X = t 
and X = n are shown in Fig. 9. Notice that the respec- 
tive exponents D n (2 — t„) = 1 and z(2 — r t ) = 1/2 are 
exact values, thus leaving only one parameter, namely k, 
to be determined by the fitting procedure. This is an ad- 
vantage of having the exact solution for the elementary 
avalanches ^From the best fit we find k = 1.14 ± 0.1 
in the given range of values of c (see caption to Fig. 9). 
It is evident from Fig. 9 that the scaling function defined 
in Eq. ( [Tl"| ) increases faster than an exponential. 

In the simulations a continuous flow may occur in the 
case of dissipative dynamics at finite lattice size L when 
the driving rate is increased. However, with increased 
system size L the behavior is different from the one in the 
case of conservative dynamics discussed above. ^From 
the numerical simulations alone we can not distinguish if 
the affected area of a continuous avalanche diverges with 
the system size L — > oo, or it remains finite for the range 
of driving rates considered here. We will also discuss 
large L limit in Sec. IV. 

We have restricted our analysis to the case where the 
degree of dissipation is such that £ <C L. For c — > 0, 
however, we have that £ — > oo, thus the role of L and £ is 
interchanged at some finite L. In the reverse limit when 
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L <C £ the behavior is expected to be similar to the case 
of conservative dynamics at finite L, studied in Sec. II. 
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FIG. 9. Finite-size scaling plot of the average duration (t) 
(top) and mass (n) (bottom) according to Eq. (pd|). Data are 
taken for values of c=0. 01,0. 0125,0.015, and 0.02 such that 
£ < L — 192 is satisfied. Note that the exponents are exact 
and the fitted value k = 1.14 within numerical error bars. 



IV. ANALYTICAL RESULTS 

We start the analytical description of the model with 
the assumption that individual avalanches, which form 
each combined avalanche, are statistically independent. 
By the definition of the model, each toppling in a sin- 
gle avalanche occurs later than those in the previous 
avalanches, so that the individual avalanches never in- 
tersect in the space-time points. Nevertheless, the next 
avalanche is sensitive to the configuration of occupation 
numbers left by the previous avalanches. In this way 
the individual avalanches are dependent on preceding 
avalanches. On the other hand, it is known from Abelian 
properties of the directed sandpile model m, that the 
recurrent state is characterized by the independent dis- 
tribution of occupation numbers zero and one at each 
site. Hence, one can expect that this property of re- 
current state provides independent distribution of single 
avalanches at least for asymptotically large systems. This 
assumption allows us to consider the process of driving 
of the directed sandpile automaton as a sequence of inde- 
pendent identically distributed (i.i.d.) events. Another 
important consequence of the statistical independence is 
that stops of combined avalanches can be considered as 
recurrent events, i.e., the probability of two successive 
stops Prob{t\,t2) at the moments t\ and t\ + ti is given 



by the product Prob(t\)Prob(t2). 

We consider the probability distribution F(t) that a 
combined avalanche starting at the moment t = stops 
for the first time at the discrete moment t. This means 
that F(t) is the probability that the stop of all preceding 
elementary avalanches occurs till the moment t. Thus, 
F(t) coincides with the probability distribution of dura- 
tion of combined avalanches, when an ensemble of events 
is considered. Along with F(t), it is useful to consider the 
function U(t) which is defined as the probability that all 
preceding single avalanches stop to the moment t regard- 
less of how many stops of combined avalanches occurred 
before t. Noting that stops of combined avalanches are 
recurrent events, we can write for F(t) and U(t) the fol- 
lowing identity 

U(t) = F(l)U(t - 1) + F(2)U(t - 2) + ... + F(t)U(Q) 

(12) 

where it is convenient to put F(0)=0 and U(0)=1. For 
the generating functions defined by 



t=0 



and 



f( S )=j2^(t) 



(13) 



(14) 



t=i 



one easily gets from Eq. ( |l2|) the known equation of the 
theory of recurrent events [EfJ 



u(s) - 1 
u(s) 



(15) 



The total probability that a combined avalanche ever 
stops is given by /(l). Therefore, the probability that 
a combined avalanche never stops, i.e., the probability of 
a continuous flow, which is given by 

F(oo) = 1-/(1), (16) 

does not vanish if u(l) in Eq. (Ill) is finite u(l) < oo . 



A. Case c = 

In the case of conservative dynamics (c = 0) the prob- 
ability distribution of durations of elementary avalanches 
is given by Eq. (|l|) . Then we can estimate the probability 
U(t) as follows. Let At = 1/r be the average time in- 
terval between addition of successive particles to the first 
row. The probability Prob(x < t) that a single avalanche 
has duration less than t is 



Prob(x < t) ~ 1 - 



(17) 
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for large t, where b is a constant of the order unity when 
L is large. Then for times t ^> At we have 



U(t) 



1 



b 



1 



(t - At) 1 / 2 



1 



(A*)!/2 

(18) 



Introducing k = t/ At we can write Eq. (|18| ) as the sum 

k ' fc 1 / 2 

n=l 



lnC/(t) = J^ln ( 1 



(nAt)V2 ; " (At)V2 



(19) 



Approximating the sum by an integral leads to 

U{t) ~ ) exp [-2ft(»Vt - VF)] (20) 

for t 3> 1 and < r < 1. For the infinite lattice there 
exists a constant C2 such that 



t(l) < C2 exp(— rVi) < oo 



(21) 



Therefore, we can find that for all finite driving rates 
r > there is a non-zero probability of continuous flow. 
In particular, the sum in Eq. (|2ll) diverges at small r as 



u(l) 



1 



(22) 



leading to the probability of stop /(l) ~ 1 — 2r 2 6 2 , which 
decreases from unity as soon as a finite driving rate is ap- 
plied. Then the probability of continuous flow Eq. (pT 
increases from zero by the same amount, i.e., 



F(oo) - r 2 ,r -> . 
For large r we expect that 

1 — F(oo) ~ exp(— r) 



(23) 



(24) 



If the size of the system is finite (L < oo) the probability 
of stops U(t) is bounded from above by a finite value 



U(t) < (rL)- b2r / 2 exp{-2brL 1 / 2 ) 



(25) 



which follows from Eq. ( p0| ) taking only the dominant 
t-behavior for t — L 3> 1. 



B. Case c > 

In the case of finite dissipation rate c > the dis- 
sipation leads to a finite characteristic length £ in the 
distribution of elementary avalanches in Eq. (|^). This 
can be easily demonstrated using mean-field arguments 
pjj in reaction-diffusion systems. Let us suppose that 



at each site of the lattice one of species A or B is liv- 
ing. These species represent two possible states of the 
original model: A corresponds to the empty site, B to 
the occupied site. Due to the external driving force, new 
particles <j> are added to the first row of the lattice at rate 
r. The propagation of particles can be described by the 
following rules 



A- 



B, B + (f>^ A + 2cf) . 



(26) 



The kinetic equations corresponding to this scheme of 
"chemical" reactions are 

h A (l) = n4l)[n B (£)-n A {£)] (27) 



h B {£)=n^e)[n A (e)-n B {£)] (28) 

M £ ) = + - 1)«bW(1 " c) (29) 

where ua{£), tib(£), and n^(£) are concentrations of 
species A, B and </>, respectively, at the £-th row. In the 
steady state we have tia = tlb — — and Eqs. (p7|)- 
( p9| ) lead to the simple conditions for the concentrations 
n A = n B = 1/2 and H 



(30) 



For c > 0, the density of particles <fi and, hence, the num- 
ber of topplings in an avalanche decay exponentially with 
the distance I from the top as n^,{£) = rexp(— £/£). This 
implies that the characteristic length of the avalanche is 
~ — ln(l — c) ~ c. Therefore the above results, in 
particular Eq. (|2^), obtained for the case of finite lat- 
tices L and c = apply also for the case c > by the 
substitution L — > L c with 

L c = min{£, L} . (31) 



C. Bounds for the busy time 

If the combined avalanches are finite (i.e., there is no 
continuous flow), we can estimate their average duration 
using the known theorem from the theory of recurrent 
events (Ref. @, Ch. XIII, Theorem 3 ). According 
to this theorem, the inverse average time of combined 
avalanches (t)^ 1 coincides with the limit of the sequence 
U(t) when t — > oo. Using the bound for U(t) given by 
Eq. H), we get 



(t) > (rLf r / 2 eM^brL 1 / 2 ) 



(32) 



The true asymptotics of (t) possibly contains an addi- 
tional prefactor L 1 / 2 like in the case for r = 0. Notice 
that we get numerically that the average duration as a 
function of r (cf. inset to Fig. 7) increases faster than 
the exponential, which agrees with Eq. (32). 
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The combination rL 1 / 2 appears as a characteris- 
tic parameter determining the duration of combined 
avalanches. Thus, for £ < L it follows from Eqs. j3^ ) 
and ( |3(]| ) that the coherence length remains constant if r 
varies with casr~ y/c, i.e. the increasing driving rate 
compensates the dissipation. 

Another interesting feature of the probability distri- 
butions at hnite driving rates is the occurrence of a 
stretch-exponential cutoffs both in dissipative and non- 
dissipative case (cf. Figs. 3,4,7). Indeed, we can see from 
Eq. @ that U(t) > F(t) for all finite t. Therefore, 
for non-dissipative case we have an exponential decay of 
combined avalanches 



V. CONCLUSIONS AND DISCUSSION 



P(t) ~ F(t) < (rt)- b2r / 2 exp(~2brVi) 



(33) 



in the thermodynamic limit L — > oo, which follows di- 



rectly from Eq. (20) for large t. 

For finite lattice sizes L or finite dissipation c > the 
function U(t) is bounded from above by a constant given 
in Eq. ( p5| ) with L replaced by L c . In this case we can 
find the origin of an exponential cut-off in the following 
way. Consider an enveloping process which corresponds 
to propagation of the front of combined avalanches. Du- 
ration of an elementary avalanche starting at tio is a sim- 
ple linear function of time, t% = t — £jo- The enveloping 
process consists of those topplings which occur at the 
maximal distance from the time axis at each moment of 
time t. The position x of the front is an one-dimensional 
random walk confined to the interval [0, £]. Starting from 
x = 0, the walk performs a step ahead with some proba- 
bility, and a step back, the length of which is a random 
variable. The combined avalanche stops if the random 
walk returns to the origin where it is trapped. The prob- 
ability of the return to the origin from an arbitrary po- 
sition x is not smaller than U (£) , where 



= K) 



-b 2 r/2 



exp(— 26r£ 



l/2\ 



(34) 



and £ ~ 1/c as above. The survival time t of the random 
walk under consideration does not exceed the period of 
successful tests in the Bernoulli scheme with the proba- 
bility of "success" 1 — U(£). The period of tests in the 
Bernoulli process has the geometrical distribution 



Prob(x = *) = (!- [7(0)*C/(0 



(35) 



Hence, the time distribution for combined avalanches is 
bounded from above by the exponential function in Eq. 
©. Using Eq. (fj), we obtain 

P(t) ~ F(t) < A(l - U(0Y ~ Aexp(-tU(0) , (36) 



where A is a constant and E/(£) is given by Eq. (34). 
Note that this expression represents an upper bound of 
the distribution of busy periods in Eq. (0). Therefore 
l/t/(£) plays the role of an effective correlation length 
at finite r, in a qualitative agreement with the numerical 
data and Eq. (|§) of previous Section. 



We have shown that a finite driving rate r is a relevant 
perturbation, which alters self-organized critical states 
in the directed sandpile automata. In the case of conser- 
vative dynamics r couples to (t)o ~ L 1 ^ 2 , thus leading 
to enhanced effects when the length scale is increased 
L — > oo. A continuous flow eventually occurs, in which 
critical long-range correlations are destroyed. On a finite 
length scale, L c = min{£, L} < oo either due to finite 
screening length £ or finite system size L, the critical 
states occur with qualitatively new correlation proper- 
ties, which is manifested in (i) a multifractal scaling of 
combined avalanches when L -C £, and (ii) occurrence 
of compensation between driving and dissipation along 
a line r (c) ~ ^ 1 / 2 ^ y/c, when I > ^. How pre- 
cisely the effective coherence length £ e / j (r) of combined 
avalanches increases with driving rate depends on details 
of the relaxation process and grain addition. In the case 
of a finite input current at each site of the system, we 
find a finite toppling rate at all scales, n^(£) ~ r/c, com- 
patible with £ e ff{r) ~~ > 00 ■ However, if grains are added 
only at the top, the correlation length increases expo- 
nentially with r in the range < r < 1, but remains 
finite presumably up to large driving rates. Here we re- 
stricted the discussion to the case r < 1, where queues of 
Dhar-Ramaswamy avalanches occur. Avalanche queuing 
for this range of driving rates is peculiar to our model, 
due to strictly local critical height rule and the directed 
transport. In the rice-pile and in the symmetric Abelian 
models fl^,[l3| a perfect queuing is prevented by the col- 
lision of avalanches, which occurs at any finite r > 0. 
Owing to the exact solution for behavior of elementary 
avalanches ||, we were able to study properties of the 
queue in detail. In particular, we find that average 
busy period is bounded from below by an exponential 
function (t) > {rVf r l' 2 exp (2brL 1 / 2 ). The avalanche 
queue can be regarded as a multifractal set, in which 
the average length is regulated by the driving rate as 
(N) = 1 + rL 1 / 2 . There are no waiting times for ele- 
mentary avalanches, therefore our model corresponds to 
the realization known in the queue theory as "infinite 
number of servers" . Hence, the average number of jobs 
q n that can be served in parallel is unlimited. It should 
be stressed that our cellular automaton represents a new 
example in the queue theory in which queuing jobs are 
distributed according to a power-law distribution with 
the exponent v < 2 and average duration of jobs is lim- 
ited by a control parameter L c = min{£, L}. We hope 
that more practical examples of this class can be found. 
We also believe that the study of the scaling properties 
of the queues, as we have done in this work, adds a new 



10 



aspect which has not been considered so far in the queue 
theory. The observed nonuniversal scaling properties of 
avalanche queues can be related to variation of the av- 
erage length of the queue with driving rate. The scaling 
exponents are found to vary approximately linearly with 
the driving rate r. A similar r-dependence was observed 
experimentally in other driven disordered systems, and 
seems to apply more generally. 

A continuous activity on the lattice, corresponding to 
a flow phase occurs for r > L^ 1 / 2 in the case of conserva- 
tive dynamics. On an infinite lattice L — > oo probability 
of continuous flow increases from zero as F(oo) ~ 2b 2 r 2 , 
whereas probability of an intermittent avalanche- like flow 
decreases from unity with the same rate /(l) ~ 1 — 2b 2 r 2 . 
For the case of dissipative dynamics c > an effective co- 
herence length increases with r for 1/c < L — > oo. Our 
results suggest that the probability /(l) remains finite 
even at high driving rates. In the limit L — > oo the com- 
pensation line extends to the point c — * 0, r — > 0. For 
the range of driving rates studied in this work we expect 
that the transport properties of grains in this model re- 
main unchanged (looked at the time scale of avalanche 
propagation) , compared to the transport at zero driving 
rate f23|| . Collision of avalanches, which occurs first at 
rates r > 1 may accelerate the grain transport, possibly 
resulting in a new scaling behavior of the distribution 
of transit times. Notice that due to local critical height 
rules and deterministic topplings the depth of the active 
zone (defined in Ref. does not change in the flow 

phase of our model. 

The analytical results in Sec. IV are derived assuming 
that elementary avalanches may be considered as inde- 
pendent events. We checked by computing numerically 
correlation function between events in a queue for finite 
L that rather weak correlations occur. The correlations 
increase with the "distance" r between avalanches as t v , 
where r\ — 0.05 ± 0.01 with statistical error bars. 
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